rm(list = ls())

## SET YOUR FILE PATH
setwd("~/Dropbox/Immigration/Analysis/Policy Diffusion/Replication Files/")

library(dplyr)
library(survminer)
library(tidyr)
library(standardize)
library(survival)
library(coxme)
library(frailtySurv)
library(frailtyEM)
library(frailtypack)
library(Hmisc)
library(haven)
library(foreign)
library(coefplot)
library(stargazer)
library(simPH)
library(coefplot)
library(xtable)
library(ggplot2)
library(texreg)
library(sp)
library(sf)
require(rgdal)
library(tmap)
library(tmaptools)
library(sp)
library(leaflet)
library(haven)
library(ggplot2)
library(shinyjs)
library(RColorBrewer)

data <- read_dta("dwrap_index")

data.africa <- subset(data, africa == 1)
data.1980 <- subset(data, year >= 1980)
data.1990 <- subset(data, year >= 1990)
data.indp5 <- subset(data, indp5 == 1)
data.indp10 <- subset(data, indp10 == 1)

set.seed(8675309)

## Set Convergence Criteria
coxph.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000, toler.inf = sqrt(.Machine$double.eps), outer.max = 1000)

coxme.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000)

## DISAGGREGATED POLICY LIBERALIZATION

m3a <- coxph(formula = Surv(accesslib1_start, accesslib1_end, accesslib1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               majorannintra_contig_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3a)

m3b <- coxph(formula = Surv(serviceslib1_start, serviceslib1_end, serviceslib1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               majorannintra_contig_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3b)

m3c <- coxph(formula = Surv(lhoodlib1_start, lhoodlib1_end, lhoodlib1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               majorannintra_contig_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3c)

m3d <- coxph(formula = Surv(movelib1_start, movelib1_end, movelib1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               majorannintra_contig_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3d)

m3e <- coxph(formula = Surv(particlib1_start, particlib1_end, particlib1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               majorannintra_contig_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3e)

m3f <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               majorannintra_contig_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3f)

disaggplot <- multiplot(m3a, m3b, m3c, m3d, m3e, m3f, secret.weapon=T, innerCI=2, coefficients = c("exc1inc2_reg_lag1"), sort = c("magnitude"), names=c("Field 1: Access \n +1 SD Liberalization", "Field 2: Services  \n +1 SD Liberalization", "Field 3: Livelihoods \n +1 SD Liberalization", "Field 4: Movement \n +1 SD Liberalization", "Field 5: Participation \n +1 SD Liberalization", "Full Index \n +1 SD Liberalization"), horizontal = F, color = c("black", "black" , "black" , "black" , "black"), zeroColor ="red", textAngle = 0, title= "", xlab = "Standardized Coefficient  \n", ylab="")
disaggplot

